// --------------------------------------------------------------------
//
//  File:        OAk8D1.rtf [Optimal Ate pairing for Jacobi quartic curves with quartic twists]
//  Date:        Oct. 10, 2018
//  Last update: Oct. 10, 2018
//  Description: Magma implementation of the optimal ate pairing on Jacobi quartic curves with quartic twists.
//               Embedding degree k = 8, CM discriminant D = 1
//
//  (c) 2018, Georgios Fotiadis & Chloe Martindale
//                 
//  Georgios Fotiadis, University of the Aegean, Greece, email: gfotiadis@aegean.gr
//  Chloe Martindale, Technische Universiteit Eindhoven, the Netherlands, email: chloemartindale@gmail.com
// --------------------------------------------------------------------
//  Polynomial Family
//  p(x) = (x^8 + x^6 + 5*x^4 + x^2 + 4*x + 4)/4
//  t(x) = x^4 + x + 2
//  r(x) = x^4 + 1 
// --------------------------------------------------------------------
//  Optimal ate pairing: a(Q,P) = [ f_{x,Q}(P) ]^((p^k - 1)/q)
// --------------------------------------------------------------------

//Miller's function point addition: h_{P,S}(Q) 
hADD:=function(P, S, d, M, N)
  U:=(P[2] + P[5]);
  V:=(S[2] + S[5]);
  R1:=P[1];
  S1:=(P[3]*S[1]);
  A:=(S1 - R1);
  A:=(A*V);
  A:=(A*U);
  U:=(S[4]*U);
  V:=(P[4]*V);
  B:=(R1*V - S1*U);
  D:=(P[1]*S[1]);
  E:=(d*D^2);
  D:=(D*(U - V));
  X3:=((R1 + S1)*(R1 - S1));
  W1:=(P[1]*P[3]);
  W2:=S[1];
  Z3:=(W1*S[2] - W2*P[2]);
  U:=(P[2]*S[2]);
  V:=P[3];
  V:=(V^2 + E);
  E:=((R1 - S1)^2);
  U3:=(X3^2);
  V3:=(Z3^2);
  Y3:=(E*(U + V) - V3);
  H:=A + D*M + B*N;
  R:=Vector(5,[X3,Y3,Z3,U3,V3]);
  return R, H;
end function;

//Miller's function point doubling: h_{P,P}(Q) 
hDBL:=function(P, d, M, N)
  U:=(P[4]^2);
  V:=(P[5]^2);
  Z3:=(V - d*U);
  E:=(P[1]*P[3]);
  D:=(2*P[4]*E);
  A:=(P[2]^2 + P[2]*P[5]);
  B:=(-P[4]*(P[2] + 2*P[5]));
  X3:=(2*E*P[2]);
  V3:=(Z3^2);
  Y3:=(2*V - Z3);
  Y3:=(2*Y3^2 - V3);
  U3:=(X3^2);
  H:=A + D*M + B*N;
  R:=Vector(5,[X3,Y3,Z3,U3,V3]);
  return R, H;
end function;

//convert Weierstrass point to extended projective Jacobi point 
ExtJacobi:=function(P)
  X:=(2*P[1])/P[2];
  Y:=(2*P[1]^3 - P[2]^2)/P[2]^2;
  Z:=(P[2])/P[2];
  U:=(X^2);
  V:=(Z^2);
  ExP:=Vector(5,[X,Y,Z,U,V]);
  return ExP;
end function;

//final exponentiation k = 8, D = 1
finalexp:=function(f, x0, p, q, Fp)
  t0:=1;
  f:=Frobenius(f, Fp, 4)/f; 
  tmp:=x0^2;
  l:=(tmp^2 + tmp + 4) div 4;
  y:=f^(l);
  t:=Frobenius(y, Fp, 3);
  t0:=t0*t;
  y:=y^x0;
  t:=Frobenius(y, Fp, 2);
  t0:=t0*t;
  y:=y^x0;
  t:=Frobenius(y, Fp, 1);
  y:=y^x0;
  t0:=t0*t*y*f;
  return t0;
end function; 

//Optimal ate pairing calculation
OptAtePairing:=function(Q, q, d, p, k, z0, M, N, Fp)
  f:=1;    
  T:=Q;
  s:=Intseq(z0,2);
  for i:=#s - 2 to 0 by -1 do
    T,h:=hDBL(T, d, M, N); 
    f:=h*f^2; 
    if s[i+1] eq 1 then
      T,h:=hADD(T, Q, d, M, N); 
      f:=h*f; 
    end if;
  end for;
  f:=finalexp(f, z0, p, q, Fp);
  return f;
end function;

//------------
//MAIN PROGRAM
//------------
t2:= 0 ; LOOP:= 100 ;
D := 1 ; k := 8 ;                 
p := 3351951982486667453837338848452726606028033606935065896469552348842908133596080151967071453287452469772970937967942438575522391344438242727636910570385409 ;
q := 115792089237333781753527356399883033864393269599672399234492847597502727192577 ;      
h := 28948022309333445438381839099970758466183387991648340884686404117030921830400 ;    
z0:= 18446744073710252032 ; 
a:= 1 ; 
b:= 0 ;

s:=k div 4; 

Fp := GF(p); 
P<x> := PolynomialRing(Fp); 
E:=EllipticCurve([Fp!(a), Fp!(b)]); 
n:=h*q; 
t:=p + 1 - n; 
f:=x^3 + a*x + b; 
d:=Fp!(-a/4); 

P<z> := PolynomialRing(Fp);
g:=IrreduciblePolynomial(Fp,s); 
Fps<u> := ExtensionField<Fp,z|g>; 
R<zz> := PolynomialRing(Fps);

for i:= 1 to 5 do
  for j:= 0 to 5 do
    fX:=i*u + j; 
    FX:=zz^(k div s) - fX;
    if IsIrreducible(FX) eq true then
      Fpk<v> := ExtensionField<Fps,zz|FX>; 
      x0:=Roots(FX, Fpk); 
      w:=x0[1][1];
      break i;
    end if;
  end for;
end for;

Et := EllipticCurve([Fps!(a*w^4),Fps!(b*w^6)]);
dt:=Fps!(d*w^4); 
ht:=#Et div q;

Ek:=EllipticCurve([Fpk!(a), Fpk!(b)]); 
Zk<xx> := PolynomialRing(Integers());
Nk := Resultant(xx^k - 1, xx^2 - t*xx + p);
hk := Nk div q^2;
dk:=Fpk!(d); 

for i:= 1 to LOOP do
//create two points PW and Pa in E(Fp) of order q (two points because of the bilinearity check)
G:=Random(E); 
nP:=h; 
P:=nP*G; 
Pa:=P;                     
nP:=Random(1, q - 1); 
PW:=nP*P; 

M:=(w^3*PW[2])/4;
N:=(w^2*PW[1])/2;

Ma:=(w^3*Pa[2])/4;
Na:=(w^2*Pa[1])/2;

//create two points QW and Qa in Et(Fp^2) of order q (two points because of the bilinearity check)
G:=Random(Et); 
nQ:= ht ; 
Q:=nQ*G;  
QW:=Q; 
xQ:=QW[1]*w^(-2); 
yQ:=QW[2]*w^(-3); 
E1:=Ek![xQ,yQ]; 
XQ:=Frobenius(xQ, Fp, 1); 
YQ:=Frobenius(yQ, Fp, 1); 
D1:=Ek![XQ,YQ]; 
QW:=D1-E1; 

xQ:=QW[1]*w^(2);
yQ:=QW[2]*w^(3);
QW:=Et![xQ,yQ]; 

nQ:=Random(1, q - 1); 
Qa:=QW; 
QW:=nQ*QW;

QJ:=ExtJacobi(QW); 
QJa:=ExtJacobi(Qa);

nPQ:=nP*nQ;

t1 := Cputime();
opt:=OptAtePairing(QJ, q, dt, p, k, z0, M, N, Fp); 
t2:=t2 + Cputime(t1); 
end for;
t2:=t2/LOOP;
printf"average time for %o pairings = %o\n", LOOP, t2;

//Bilinearity check
OptAtePairing(QJ, q, dt, p, k, z0, M, N, Fp) eq OptAtePairing(QJa, q, dt, p, k, z0, Ma, Na, Fp)^(nPQ);